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ABSTRACT 

The development and application of benchmark examples for the assessment of quasi- 
static delamination propagation capabilities was demonstrated for ANSYS® and 
Abaqus/Standard®. The examples selected were based on finite element models of 
Double Cantilever Beam (DCB) and Mixed-Mode Bending (MMB) specimens. First, 
quasi-static benchmark results were created based on an approach developed 
previously. Second, the delamination was allowed to propagate under quasi-static 
loading from its initial location using the automated procedure implemented in 
ANSYS® and Abaqus/Standard®. Input control parameters were varied to study the 
effect on the computed delamination propagation. Overall, the benchmarking 
procedure proved valuable by highlighting the issues associated with choosing the 
appropriate input parameters for the VCCT implementations in ANSYS® and 
Abaqus/Standard®. However, further assessment for mixed-mode delamination 
fatigue onset and growth is required. Additionally studies should include the 
assessment of the propagation capabilities in more complex specimens and on a 
structural level. 


INTRODUCTION 

Over the past two decades, the use of fracture mechanics has become common 
practice to characterize the onset and growth of delaminations. In order to predict 
delamination onset or growth, the calculated strain energy release rate components 
are compared to interlaminar fracture toughness properties measured over a range 
from pure mode I loading to pure mode II loading. 
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The virtual crack closure technique (VCCT) is widely used for computing 
energy release rates based on results from continuum (2D) and solid (3D) finite 
element (FE) analyses and to supply the mode separation required when using the 
mixed-mode fracture criterion [1, 2], The virtual crack closure technique was 
recently implemented into several commercial finite element codes such as 
Abaqus/Standard®, MD Nastran™, Marc™ and ANSYS®. As new methods for 
analyzing composite delamination are incorporated into finite element codes, the 
need for comparison and benchmarking becomes important since each code 
requires specific input parameters unique to its implementation. These parameters 
are unique to the numerical approach chosen and do not reflect real physical 
differences in delamination behavior. 

An approach for assessing the mode I, and mixed-mode I and II, delamination 
propagation capabilities in commercial finite element codes under static loading was 
recently presented and demonstrated for the VCCT implementation in 
Abaqus/Standard® [3-5] as well as MD Nastran™ and Marc™ [6] and ANSYS® [7]. 
First, benchmark results were created manually for finite element models of the mode 
I Double Cantilever Beam (DCB), the mode II End Notched Flexure (ENF) as well as 
the mixed-mode Eli Single Leg Bending (SLB) and Mixed-Mode Bending (MMB) 
specimens. Second, the delamination was allowed to propagate under quasi-static 
loading from its initial location using the automated procedure implemented in the 
finite element software. The approach was then extended to allow the assessment of 
the delamination fatigue growth prediction capabilities in commercial finite element 
codes [4, 8]. As for the static case, benchmark results were created manually first for 
the mode I Double Cantilever Beam (DCB) and the mode II End Notched Flexure 
(ENF) specimen. Second, the delamination was allowed to grow under cyclic loading 
in a finite element model of a commercial code. For all cases, input control parameters 
were varied to study the effect on the computed delamination propagation and growth. 
The benchmarking procedure proved valuable by highlighting the issues associated 
with choosing the input parameters of the particular implementation. Consequently, 
the benchmark enabled the selection of the appropriate input parameters that yielded 
good agreement between the results obtained from the growth analysis and the 
benchmark results. Once the parameters have been identified, they may then be used 
with greater confidence to model delamination growth for more complex 
configurations. 

In this paper, the benchmark analysis from previous studies is summarized using 
the Double Cantilever Beam (DCB) and the Mixed-Mode Bending (MMB) specimens 
as examples [3, 5, 7]. Special attention is given first to the creation of the benchmark 
case, which is independent of the analysis software. The benchmark case allows the 
assessment of the automated delamination propagation prediction capabilities in 
commercial finite element codes based on the virtual crack closure technique (VCCT). 
Examples of automated propagation analyses are then shown for two different finite 
element codes (Abaqus/Standard® and ANSYS®) and the selection of the required 
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code specific input parameters is discussed. Examples using MD Nastran and Marc 
are discussed in reference 6. 


SPECIMEN CONFIGURATIONS SELECTED AS BENCHMARK CASES 


For the current summary, two simple specimens used for fracture toughness 
testing were selected. These specimens had been used previously to develop the 
current approach for assessment of the quasi-static delamination propagation 
simulation capabilities in commercial finite element codes [3-7]. The Double 
Cantilever Beam (DCB) specimen, as shown in Figure 1, was chosen since it only 
exhibits the mode I opening fracture mode. Besides the previously developed 
benchmark case of the DCB specimen [3], an additional example was used for which 
experimental data also exist for comparison. This example was published by NAFEMS 
- an independent not-for-profit body with the sole aim of promoting the effective use 
of engineering simulation methods - as part of their benchmark cases [9]. 

One example that exhibits mixed-mode I/II fracture was also selected: The Mixed- 
Mode Bending (MMB) specimen, as shown in Figure 2, was studied for 80% mode II 
[5,7]. The material and overall specimen dimensions including initial crack length, ao, 
are shown in the respective figures. All configurations had layups of [0]24. The 
material properties are given in references 3, 5 and 7. A methodology for delamination 
propagation, onset and growth based on fracture mechanics [10] was used to create the 
benchmark examples. 
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Krueger example [3] 

8 25.0 mm 

2h 3.0 mm 
21 150.0 mm 
ag 30.0 mm 
material and layup: 
T300/1076 [0] 24 

NAFEMS example [9] 
8 30.0 mm 

2h 3.0 mm 
2L 150.0 mm 
a 0 30.0 mm 

material and layup: 
T800/924 [0] 24 


Figure 1. Double Cantilever Beam specimen (DCB): Dimensions and Layup [3,7]. 
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b 25.4 mm 
2h 4.5 mm 
2L 100.8 mm 
2L* 152.4 mm 
a 0 25.4 mm 

material and layup: 
IM7/8552 [0] 24 


Figure 2. Mixed-Mode Bending specimen (MMB): Dimensions and Layup [5,7]. 


METHODOLOGY BASED ON FRACTURE MECHANICS 


A quasi-static mixed-mode fracture criterion is discussed first, since the 
parameters are required input to the VCCT implementation in Abaqus/Standard® and 
ANSYS®. The input details of the specific codes are discussed in references 3 and 7. 
The mixed-mode fracture criterion for a material is determined by plotting the 
interlaminar fracture toughness, G c , versus the mixed-mode ratio, Gu/Gt, as shown in 
Figure 3 for a typical carbon/epoxy material. The fracture criterion is generated 
experimentally using pure Mode I {Gn/Gr=0) Double Cantilever Beam (DCB) tests 
(as shown in Figure 1), pure Mode II (G///Gt= 1) End-Notched Flexure (ENF) tests, 
and Mixed Mode Bending (MMB) tests (as shown in Figure 2), of varying ratios of 
Gii/Gt. The mean values (filled blue circles) are shown in Figure 3. A 2D fracture 
criterion was suggested by Benzeggah and Kenane [11] using a simple mathematical 
relationship between G c and Gu/Gt 

G c = Gj c + (G IIC — G lc ) (1) 

In this expression, typically called the B-K criterion, G/ c and G// c are the 
experimentally determined fracture toughness data for mode I and II as shown in 
Figure 3. The exponent rj was determined by a curve fit. The parameters G/ c , Gu c 
and t] are required input to perform a VCCT analysis in Abaqus/Standard® [12] and 
ANSYS® [13, 14]. Example input is provided in references 3, 5 and 7. 

During an automated propagation analysis, the total strain energy release rate, 
Gt, and the mixed-mode ratio Gu/Gt are computed using VCCT. The failure index, 
Gt/G c , is calculated by correlating the computed total energy release rate, Gt, with 
the mixed-mode fracture toughness, G c , of the graphite/epoxy material. As shown in 
equation (1) the mixed-mode fracture toughness, G c , is a function of the mixed- 
mode ratio Gu/Gt ( see also Figure 3). An example how to obtain G c for Gii/Gt=0.S 
has been included in Figure 3 (purple arrows). It is assumed that the delamination 
propagates when the failure index, Gt/G c , reaches unity. 



Figure 3. Mixed mode fracture criterion for a typical carbon/epoxy material. 


FINITE ELEMENT MODELING 


Model description 

Finite element models of the specimens are shown in Figures 4 and 5. For two- 
dimensional (2D) analyses the specimens were modeled with plane strain elements 
such as CPE41 in Abaqus/Standard® [12] and PLANE182 in ANSYS® [13-14]. 
Three-dimensional (3D) analyses were only performed using Abaqus/Standard® and 
solid brick elements C3D8I were used to model the specimens. Along the length, all 
models were divided into different sections with different mesh refinement as 
shown in Figures 4 and 5. A typical finite element model of a DCB specimen is 
shown in Figure 4. The DCB specimen was modeled with six elements through the 
specimen thickness (2h). The resulting element length at the delamination tip was 
Aa= 0.5 mm. 



(a). Deformed model of a DCB specimen ( plane strain CPE4 elements in 
Abaqus/Standard® and PLANE! 82 elements in ANSYS®) [7,8]. 





(b). Three-dimensional model (solid C3D8I elements in Abaqus/Standard® )[3 ] . 


Figure 4. Double Cantilever Beam specimen (DCB). 
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Figure 5. Two-dimensional model ( plane strain CPE4I elements in A baqus/Standanf ’ ) 

of a MMB specimen [5,7]. 


A two-dimensional finite element model of a MMB specimen with boundary 
conditions is shown in Figure 5. Along the length, all models were divided into 
different sections with different mesh refinement as discussed earlier. The MMB 
specimen was also modeled with six elements through the specimen thickness. The 
resulting element length at the delamination tip was Aa=0.5 mm. The load 
apparatus was modeled explicitly using rigid beam elements (R2D2) as shown in 
Figure 5. Multi-point constraints were used to connect the rigid elements with the 
solid model of the specimen and enforce the appropriate boundary conditions. The 
mixed-mode ratio Gu/Gt is controlled by the length, c, of the loading arm. 
Configurations were developed that yielded mode ratios of 20% mode II ( Gu/Gt 
=0.2), 50% mode II ( G U /G T =0.5), and 80% mode II (G„/G T =0.8) [5,7], Here, only 
the model for 80% mode II is shown for which results will be presented later. 

For all models, the plane of delamination was modeled as a discrete 
discontinuity in the center of the specimen. The models were created as separate 
meshes for the upper and lower part of the specimens with identical nodal point 
coordinates in the plane of delamination and contact was used to define the intact 
section of the specimen. Details of the modeling are discussed in references 3, 5 and 
7. 


Quasi-static delamination propagation analysis 

For the automated delamination propagation analysis, the VCCT 
implementation in Abaqus/Standard® 6.10 [12] and in ANSYS® 13.0 beta and 14.0 
beta were used [13, 14], It is implied that the energy release rate at the crack tip is 
calculated at the end of a converged increment. Once the energy release rate 
exceeds the critical strain energy release rate (including the user-specified mixed- 
mode criteria as shown in Figure 3), the node at the crack tip is released in the 
following increment, which allows the crack to propagate. For automated 
propagation analysis, it was assumed that the computed behavior should closely 
match the benchmark results created below. For all analyses, the elastic constants 
and the input to define the fracture criterion were kept constant. Only code specific, 


implementation dependent parameters were varied to improve convergence and 
achieve better results. 

In addition to the mixed-mode fracture criterion, Abaqus/Standard® requires 
extra input for the propagation analysis. If a user specified release tolerance is 
exceeded in an increment ( G-G c )/G c > release tolerance ), a cutback operation is 
performed which reduces the time increment. In the new smaller increment, the 
strain energy release rates are recalculated and compared to the user specified 
release tolerance. The cutback reduces the degree of overshoot and improves the 
accuracy of the local solution [12]. A release tolerance of 0.2 is suggested in the 
handbook [12]. Further, to help overcome convergence issues during the 
propagation analysis, Abaqus/Standard® provides: 

• contact stabilization which is applied across only selected contact pairs and used 
to control the motion of two contact pairs while they approach each other in 
multi-body contact. The damping is applied when bonded contact pairs debond 
and move away from each other [12] 

• automatic or static stabilization which is applied to the motion of the entire 
model and is commonly used in models that exhibit statically unstable behavior 
such as buckling [12] 

• viscous regularization, which is applied only to nodes on contact pairs that have 
just debonded. The viscous regularization damping causes the tangent stiffness 
matrix of the softening material to be positive for sufficiently small time 
increments. Viscous regularization damping in Abaqus/Standard® is similar to 
the viscous regularization damping provided for cohesive elements and the 
concrete material model in Abaqus/Standard® [12]. Further details about the 
required input parameters are discussed in references 5 and 7. 

The VCCT implementation in ANSYS® does not require the user to specify a 
release tolerance. Also stabilization criteria to help overcome convergence issues 
were not available in the ANSYS® releases studied. The following specific input 
parameters could be altered by the user in order to achieve better results during an 
automated delamination propagation analysis: 

• The initial time step size and maximum step size for crack growth ( cgrow ) were 
varied. The parameters are called dtime and dtmax. 

• The automated time stepping ( autots ) was also studied in detail. The input 
parameters, which were varied are also called dtime and dtmax 3 . Further details 
about the required input parameters are discussed in reference 9. 


PROCEDURE FOR DEVELOPING QUASI-STATIC BENCHMARK CASES 

Based on the approach developed earlier [3], quasi-static benchmark results can 
be created for any analysis software used. The procedure is outlined using the DCB 
specimen with a unidirectional layup (as shown in Figure 1) as an example. This 
specimen configuration was chosen as a benchmark case for mode I, since it is 


3 

Note that the input parameters for both crack growth time step (cgrow) and automated time stepping 
(autots) are called dtime and dtmax. For more clarity dtimea and dtmaxa were used for the parameters 
for automated time stepping (autots). 


simple and a number of numerical studies had been performed previously to 
evaluate the critical strain energy release rates. To avoid unnecessary 
complications, experimental anomalies such as fiber bridging were not addressed 
[3]. 

• First, finite element models of the specimen with different delamination lengths, 
cio, have to be created. For the current example two-dimensional finite element 
models simulating the DCB specimens with 17 different delamination lengths ao 
were created (30.5 mm<ao<69.5 mm). 

• For each delamination length, ao, modeled, the load, P, and opening 
displacement, 6/2, at the load point were plotted as shown in Figure 6 (dashed 
colored and solid grey lines) 

• For each delamination length, ao, modeled, the mixed-mode strain energy release 
rate components were computed using VCCT for applied displacements 
6/2 =1.5 mm (for ao<45.4 mm), and 6/2 =3 mm (for 45.4 mm< ao<65.4 mm) 
and 6/2 =5 mm (for 65.4 mm< ao< 69.5 mm). The total energy release rate is a 
function of the delamination length, ao, and the applied opening displacement, 
6/2, as indicated in Figure 6, thus Gr= Gp (ao, 6/2). For the simple case of the 
mode I DCB specimen shown, the computed values were predominantly mode I 
as expected. 

• For each delamination length, ao, modeled, a failure index, Gj/G c , is calculated 
by correlating the computed total energy release rate, Gp, with the mixed-mode 
fracture toughness, G c , of the material shown in Figure 3. For the simple case of 
the DCB specimen shown, the failure index is simply calculated as G//G/ c . It is 
assumed that the delamination propagates when the failure index reaches unity. 

• Therefore, the critical load, P cnh and critical opening displacement, 6 cr i t /2 can be 
calculated - for each delamination length, ao, modeled - based on the relationship 
between load, P, and the energy release rate, G [15], 


P 2 dC P 

G ~Y'~dA 


( 2 ) 


In equation (2), Cp is the compliance of the specimen, and DA is the increase 
in surface area corresponding to an incremental increase in load or displacement 
at fracture. The critical load, P cnh and critical opening displacement, & nf /2 can 
be calculated for each delamination length, ao, modeled 
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and the results can be included in the load/displacement plots as shown in 
Figure 7 (solid red circles). 

• By fitting a curve through these critical load/displacement results (solid red 
circles), a benchmark solution (thick solid red line) can be created as shown in 
Figure 7. 


During the automated propagation analysis, the computed load/displacements results 
are expected to follow the benchmark solution. 
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applied opening displacement 8 / 2 , mm 

Figure 6. Computed load-displacement behavior of a DCB specimen for different 

delamination lengths ao. 
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Figure 7. Calculated critical behavior and resulting benchmark case for a DCB 

specimen. 


QUASI-STATIC ANALYSIS BENCHMARKING 


Quasi-Static Benchmark Case for Mode I 

Automated delamination propagation analysis using Abaqus/Standarcf 

A set of example results using Abaqus/Standard® are shown in Figure 8 where the 
computed resultant force (load P) at the tip of the DCB specimen is plotted versus the 
applied crack tip opening (6/2) for different input parameters [3]. To overcome the 
convergence problems, the methods implemented in Abaqus/Standard® were used 
independently to study the effects [12]. First, global stabilization was used in the 
analysis. For a stabilization factor of 2x1 O' 5 , the stiffness changed to almost infinity 
once the critical load was reached causing the load to increase sharply (thick solid 
red). The load increased until a point was reached where the delamination propagation 
started and the load gradually decreased following a saw tooth curve with local rising 
and declining segments. The gradual load decrease followed the same trend as the 
benchmark curve (solid grey line) but is shifted toward higher loads. For a 
stabilization factor of 2x1 O' 6 (solid green line), the same saw tooth pattern was 
observed but the average curve was in good agreement with the benchmark result. 
Second, contact stabilization was used in the analysis. For all combinations of 
stabilization factors and release tolerances, a saw tooth pattern was observed, where 
the peak values were in good agreement with the benchmark result (solid blue line). 
Third, the use of viscous regularization to overcome convergence problems was 
studied. Convergence could not be achieved over a wide range of viscosity 
coefficients when a default release tolerance value of 0.2 was used. 
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Figure 8. Computed critical load-displacement behavior for a DCB specimen 
obtained from 3D solid models ( C3D8I). 



Subsequently, the release tolerance value was increased. For all combinations of 
the viscosity coefficient and release tolerance, where convergence was achieved, a 
saw tooth pattern was obtained, where the peak values were in good agreement with 
the benchmark result (thin solid black line). 

The examples shown highlight the importance of benchmarking to identify 
critical analysis input parameters. In summary, good agreement between the load- 
displacement relationship obtained from the propagation analysis results and the 
benchmark results could be achieved by selecting the appropriate input parameters. 
However, selecting the appropriate VCCT input parameters such as release 
tolerance, global or contact stabilization and viscous regularization, was not 
straightforward and often required an iterative procedure. The default setting for 
global stabilization yielded unsatisfactory results although the analysis converged. 
Detailed results for Abaqus/Standard® are discussed in reference 3. 


Automated delamination propagation analysis using A/VSYS® 

The automated propagation analysis in ANSYS® was performed in one step. 
Starting from an initial delamination length, ao=30.5 mm, the delamination was 
allowed to propagate based on the algorithms implemented into ANSYS® [7]. A 
total crack opening displacement 612=2.0 mm was applied to each cantilever arm. 
Results from the first set of propagation analyses are shown in Figure 9. 
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Figure 9. Computed critical load-displacement behavior for a DCB specimen 
obtained from 2D planar models (PLANE182). 


For the first propagation analysis (dashed red line), the stiffness of the specimen 
remained unchanged once the critical point was reached and load and displacement 
kept increasing. Later, the stiffness decreased as the delamination propagated, the 
load started to drop and the computed load/displacement path converged to the 
benchmark result (solid grey circles and solid grey line). In order to minimize the 
undesirable overshoot and to closely capture the critical point, the parameters 
(dtime and dtmax) had to be specified to control the time step for crack growth 
(cgrow). For an initial time step, dtime= 0.001, and maximum allowable time step, 
dtmax= 0.01 - so that the analysis could adjust the parameter as needed - the result 
(solid blue line) was in good agreement with the benchmark result. Improved 
results (solid red line) could be obtained when the maximum allowable time step 
was limited to dtmax= 0.001 

Additionally, the effect on the computed results caused by the input parameters 
(dtimea and dtmaxa ) for automated time stepping ( autots ) was studied. The effect 
was studied in more detail while - based on the study from above - the input to 
control the time step for crack growth was kept constant (dtime=dtmax= 10' 6 ). For 
automated time steps dtimea=dtmaxa= 0.1 (dtmina= 10' 6 ), the results (solid green 
line) overshot the critical point as shown in Figure 9. Once delamination 
propagation started, the load instantly dropped and the analysis continued to follow 
the benchmark (solid grey circles and solid grey line). For dtimea=dtmaxa= 0.01, 
the overshoot was less pronounced (solid purple line). It appears that crack 
propagation starts at the time step after which the critical point has been exceeded 
for the first time. As shown in the plots, this time step depends on the step size 
(green dots and purple diamonds). There is no apparent automated cut back 
procedure that forces a reanalysis of the last step with a smaller time step, which 
would allow the analysis to zero in on the critical point. The implementation of a 
user defined release tolerance, which triggers a cut back might alleviate this 
observed problem. 

In summary, input parameters for automated time stepping ( autots ) need to be 
chosen by the user such that the time steps are sufficiently small to capture the 
critical point for propagation onset (first node release) correctly. The input 
parameters, which control the time step for crack growth however, appear to have 
an effect only on the quality of results during propagation. The results suggest that 
small time steps are required to obtain accurate results. However, such small steps 
cause an increase in computation time as discussed later. Long computation times, 
however, may be avoided by splitting the analysis in two parts. During the first part, 
the analysis remains subcritical and the delamination does not propagate so that 
large time steps can be used for automated time stepping. The first part of the 
analysis is set to end just before the analysis reaches the critical point. During the 
second part of the analysis, the time steps for automated time stepping are reduced 
to accurately capture the critical point. Additionally, the time steps for crack growth 
are selected to be sufficiently small to assure a proper propagation analysis. 

NAFEMS DCB benchmark case for mode I 

To highlight the independence of the benchmark procedure from the analysis 
code used, another DCB specimen with slightly different dimensions and material 
properties (shown in Figure 1) was selected. For this configuration experimental as 


well as analytical results were available. This example was published by NAFEMS 
as part of their benchmark examples [8]. The experiments exhibited only a 
negligible amount of fiber bridging so that comparison with the analysis appeared 
justified [8]. Following the same procedure outlined above, two benchmark results 
were created for ANSYS® and Abaqus/Standard® using the same model shown in 
Figure 4a. As shown in the plots of Figure 10, the computed benchmark results for 
ANSYS® (open red circles and solid red line) compared well with the benchmark 
results for Abaqus/Standard® (solid black diamond and solid black line) [3] as well 
as and the experimental data (solid blue squares) published in reference 8. The 
benchmark results also agree with the analytical results (thick dashed green line). 


Quasi-Static Benchmark Cases for Mixed-Mode I/II 

Development of benchmark cases for 80% mode II based on the MMB specimen 

A benchmark case for 80% mode II ( Gu/Gt =0.8) was created based on the 
approach discussed above for mode I. Two-dimensional finite element models 
simulating MMB specimens with 15 different delamination lengths ao were created 
(25.4 mm< ao< 70.6 mm) to study the case of Gu/Gt =0.8. An example of a finite 
element model is shown in Figure 5. For each delamination length modeled, the 
load, P, and displacement, w, were monitored as shown in Figure 11 (dashed 
colored and solid grey lines). Using VCCT, the total energy release rate, Gt, and the 
mixed-mode ratio Gu/Gt were computed at the end of the analysis as shown in 
Figure 11. 
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Figure 10. Benchmark case for NAFEMS DCB specimen. 
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Figure 1 1 . Calculated critical load-displacement behavior for a MMB specimen 
(80% mode II) obtained from 2D planar models ( CPE41 ). 

The critical mixed-mode fracture toughness, G c for each computed mixed-mode 
ratio ( Gh/Gt —0.8) was obtained from the curve fit of the material data (solid red 
curve) as illustrated in Figure 3 (purple arrows). The critical load, P cr it, and critical 
displacement, w cri( , were calculated for each delamination length modeled, using 
equation (3), and the results were included in the load/displacement plots as shown 
in Figure 1 1 (solid red circles). 

These critical load/displacement results indicated that, with increasing 
delamination length, less load is required to extend the delamination. For the first 
five delamination lengths, ao, plotted in Figure 11, the values of the critical 
displacements also decreased at the same time. This means that the MMB specimen 
exhibits unstable delamination propagation under load control as well as 
displacement control in this region for Gu/Gt =0.8. The remaining critical 
load/displacement results indicated stable propagation. From these critical 
load/displacement results (solid grey circles and dashed line), two benchmark 
solutions can be created as shown in Figure 12. 

During the analysis, either prescribed displacements, w, or nodal point loads, P, 
are applied. For the case of prescribed displacements, w, (dashed blue line), the 
applied displacement must be held constant over several increments once the 
critical point ( P cr it , w cr it ) is reached, and the delamination front is advanced during 
these increments. Once the critical path (dashed grey line) is reached, the applied 
displacement is increased again incrementally. For the case of applied nodal point 
loads (dashed red line), the applied load must be held constant while the 
delamination front is advanced during these increments. Once the critical path 
(dashed grey line) is reached, the applied load is increased again incrementally. It is 
assumed that the load/displacement relationship computed during automatic 
propagation should closely match the benchmark case. Benchmark cases for 20% 
and 50% mode II are discussed in detail in reference 5. 
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Figure 12. Calculated critical behavior and resulting benchmark cases for applied 
load, P, and displacement, w, for a MMB specimen ( 80% mode II). 


Computed delamination propagation for applied displacement 

The propagation analysis was performed in two steps using the model shown in 
Figure 5 starting from an initial delamination length, ao=25A mm. In the first step, a 
displacement of w=1.5 mm was applied which nearly equaled the critical 
displacement, vv c/ «=1 .65 mm, determined earlier. In the second step, the applied 
displacement was increased to w-8.0 mm. For this second step, automatic 
incrementation was used in Abaqus/Standard® and a small increment size (0.5% of 
the total step) was chosen at the beginning of the step. To reduce the risk of 
numerical instability and early ter mi nation of the analysis, a minimum allowed 
increment size (10 18 of the total step) was also chosen. The analysis was limited to 
5000 increments. 

Initially, analyses were performed using two-dimensional planar models (shown 
in Figure 5) without stabilization or viscous regularization. Using the default value 
reltol= 0.2, the load and displacement exceeded the critical point and the analysis 
terminated (thick solid green line) as shown in Figure 13 where the computed 
resultant force (load P) is plotted versus the applied displacement w. By increasing 
the release tolerance to reltol= 0.5, - as suggested by the Abaqus/Standard® error in 
the message (.msg) file - it was possible to extend the analysis without an error 
message until the load drop occurred and delamination propagation started as 
shown in Figure 13, (solid light blue line). However, the load drop occurred past the 
critical point and the analysis terminated soon afterwards due to convergence 
problems. The release tolerance was not increased any further as suggested in the 
Abaqus/Standard® error in the message (.msg) file. Previous analysis had shown 
that by increasing the release tolerance (reltol> 0.2), termination of the analysis 
could be avoided, however, the results had not been in good agreement with the 


benchmark [3,4]. Therefore, it was decided to introduce additional stabilization to 
obtain better agreement with the benchmark case. 

The results computed when contact stabilization (cs) was added are plotted in 
Figure 14. A small stabilization factor (c.s-1 x 1 0" 6 ) was used for all cases since it 
had yielded good results in previous analyses [3, 4], Initially, the release tolerance 
value was set at the default value (reltol= 0.2) (dashed green line). For this 
parameter combination, the load increased past the critical point and delamination 
propagation started while the load dropped (dashed green line). The 
load/displacement path then ran parallel to the constant displacement branch of the 
benchmark result and later followed the stable propagation branch of the 
benchmark result. 

To reduce the observed overshoot, the release tolerance was reduced. For a 
stabilization factor of cs= 1 x 1 0‘ 6 and release tolerance reltol=0A (solid blue line), 
the overshoot was reduced, and the computed load/displacement path then ran 
closer to the constant displacement branch of the benchmark result. However, just 
before reaching the minimum load, the displacement increased slightly followed by 
a small load drop. The results then followed the stable propagation branch of the 
benchmark result. Further reducing the release tolerance (reltol= 0.01) yielded 
results that were in excellent agreement with the benchmark (dash-dotted red line). 
Due to the fine mesh, only a small saw-tooth pattern was observed along the stable 
path. In particular, good agreement could be achieved by selecting input parameters 
that had previously been determined during analyses of mode I DCB and mode II 
ENF specimens [3, 4], 
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Figure 13. Computed critical load-displacement behavior for MMB specimen (80% 
mode II) obtained from 2D planar models ( CPE4I) with different release tolerance 

settings. 
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Figure 14. Computed critical load-displacement behavior for MMB specimen (80% 
mode II) obtained from 2D planar models ( CPE4I) with added contact stabilization. 


Computed delamination propagation for applied quasi- static load 

The propagation analysis was performed in two steps using the model shown in 
Figure 5 starting from an initial delamination length, ao=25A mm. In the first step, a 
load P=130 N was applied which equaled nearly the critical load, P c m=75 1 N, 
determined earlier. In the second step, the total load was increased (/ ) =1000 N). As 
before, automatic incrementation was used with a small increment size at the 
beginning (0.5% of the total step) and a very small minimum allowed increment 
(10‘ 18 of the total step) to reduce the risk of numerical instability and early 
termination of the analysis. The analysis was limited to 5000 increments. 

The same steps discussed in the section on applied displacement were followed. 
Initially, analyses were performed using two-dimensional planar models without 
stabilization or viscous regularization. Using the default value (reltol= 0.2), the load 
increase stopped after reaching the critical point, but the analysis terminated 
immediately (thick solid green line) due to convergence problems, as shown in 
Figure 15. The release tolerance was not increased as suggested by the 
Abaqus/Standard® error in the message (.msg) file. Previous analysis had shown 
that by increasing the release tolerance (reltol> 0.2), termination of the analysis 
could be avoided, however, the results had not been in good agreement with the 
benchmark [3,4]. Therefore, it was decided to introduce additional stabilization to 
obtain better agreement with the benchmark case. 

The results computed when contact stabilization (cs) was added are plotted in 
Figure 16. A small stabilization factor (c.s-1 x 1 0" 6 ) was used for all cases since it 
had yielded good results in previous analyses [3, 4], Initially, the release tolerance 
value was set at the default value (reltol= 0.2) (dashed green line). For this 
parameter combination, the load increased up to the critical point and delamination 


propagation started while the load remained constant (dashed green line). Also for 
the stable propagation path, the result was in good agreement with the benchmark 
result (solid grey line). Then, the release tolerance was reduced to reltol= 0.1 (solid 
blue line) and further to reltol= 0.01 (dash-dotted red line). 
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Figure 15. Computed critical load-displacement behavior for MMB specimen 
( 80% mode II) obtained from 2D planar models ( CPE4I) with different release 
tolerance settings. 
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Figure 16. Computed critical load-displacement behavior for MMB specimen (80% 
mode II) obtained from 2D planar models ( CPE4I) with added contact stabilization. 


For all cases, the results were in good agreement with the benchmark results 
over the entire load/displacement range. In particular, good agreement could be 
achieved by selecting input parameters that had previously been determined during 
analyses of mode I DCB and mode II ENF specimens [3, 4]. 


SUMMARY AND CONCLUSIONS 

The development and application of benchmark examples for the assessment of 
quasi-static delamination propagation capabilities was demonstrated for ANSYS® 
and Abaqus/Standard®. The examples selected were based on finite element models 
of Double Cantilever Beam (DCB), and Mixed-Mode Bending (MMB) specimens. 

First, quasi-static benchmark results were created based on the approach 
developed in reference 3. Second, the delamination was allowed to propagate under 
quasi-static loading from its initial location using the automated procedure 
implemented in ANSYS® and Abaqus/Standard®. Input control parameters were 
varied to study the effect on the computed delamination propagation. 

The results showed the following: 

• The benchmarking procedure is independent of the analysis software. 

• The benchmarking procedure was able to recreate the DCB benchmark case 
suggested by NAFEMS. 

• The benchmarking procedure was capable of identifying the critical input 
parameters for a particular implementation and highlighting the differences in 
implementation. 

• These input parameters are unique to the software studied and do not reflect real 
physical differences in delamination behavior. 

• In general, good agreement between the results obtained from the propagation 
analysis and the benchmark results could be achieved by selecting the 
appropriate input parameters. However, selecting the appropriate input required 
an iterative procedure. 

• For Abaqus/Standard® in particular, the results for automated delamination 
propagation analysis under quasi-static loading showed the following: 

o Using the default release tolerance (reltol= 0.2) or increasing the value, as 
suggested in the user's manual, may help to overcome convergence 
problems, however, it leads to an undesired overshoot of the computed 
result compared to the benchmark. 

o A combination of release tolerance and contact stabilization is required to 
obtain more accurate results. 

o A gradual reduction of the release tolerance and contact stabilization over 
several analyses is suggested. 

o Good agreement between analysis results and the benchmarks could be 
achieved for release tolerance values (reltol<0A) in combination with 
contact stabilization (cs-lxlO" 6 ). 

• For ANSYS®, in particular, the results for automated delamination propagation 
analysis under static loading showed the following: 

o Automated analysis generally resulted in a converged solution, however, it 
led to undesired overshoot of the result. 


o User input to control the time step for crack growth ( cgrow ) is required to 
obtain more accurate results, in particular for the propagation phase of the 
analysis. Repeatedly consistent results could be obtained when the initial 
time step, chime, and the maximum allowable time step, dtmax, were 
equivalent and were equal to or smaller than 10' 4 (dtime=dtmax<\ 0‘ 4 ). 
o User input to control the automated time stepping ( autots ) is also required, 
in particular to capture the critical point of propagation onset or first node 
release. Repeatedly consistent results could be obtained when the initial 
time step, dtimea, and the maximum allowable time step, dtmaxa, were 
equivalent and were equal to or smaller than 10' 3 (dtimea=dtmaxa< 10' 3 ). 

Overall, the benchmarking procedure proved valuable by highlighting the issues 
associated with choosing the appropriate input parameters for the VCCT 
implementations in ANSYS® and Abaqus/Standard®. However, further assessment for 
mixed-mode delamination fatigue onset and growth is required. Additionally studies 
should include the assessment of the propagation capabilities in more complex 
specimens and on a structural level. 
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